cd "INSERT PATH HERE"



************************************
************************************
************************************
************************************
************************************
*Code to reproduce Figure 1
************************************
************************************
************************************
************************************
************************************
*To recreate Figure 1 load .dta file below and run code

use "County level reopening data.dta", clear

twoway scatter fullyopen ptrump [w=countypop], msymbol(Oh) || lfitci fullyopen ptrump
gr_edit yaxis1.reset_rule 0 100 25 , tickset(major) ruletype(range) 
gr_edit yaxis1.style.editstyle majorstyle(tickangle(horizontal)) editcopy
gr_edit legend.draw_view.setstyle, style(no)
gr_edit yaxis1.title.text.Arrpush % Attending In-Person Classes
gr_edit yaxis1.title.style.editstyle size(medium) editcopy
gr_edit xaxis1.title.text = {}
gr_edit xaxis1.title.text.Arrpush Trump Vote Share (2016)

graph export "Figure 1", replace as(png)






************************************
************************************
************************************
************************************
************************************
*Code to reproduce Table 1
************************************
************************************
************************************
************************************
************************************

use "Main analysis file.dta", clear


*Baseline Nationwide estimation (using state and locale FEs)
oprobit openedup i.st i.loc2 ptrump distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k, cluster(st) 
estimates store col1

*Baseline Nationwide estimation but swap to state-by-locale FEs 
oprobit openedup i.st##i.loc2 ptrump distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k, cluster(st)
estimates store col2

*Baseline Nationwide estimation but swap to Commuter Zone (CZ) FEs
oprobit openedup i.st i.loc2 i.czi ptrump distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k, cluster(st)
estimates store col3

*Baseline Nationwide stimation but use collective bargaining dummy instead of district size
oprobit openedup i.st i.loc2 ptrump cbmoe distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k, cluster(st)
estimates store col4

*Washington state (using a better measure of union political strength)
oprobit openedup i.loc2 ptrump distsize lninc50avgall logpps pwhite pac2018 meancasesper2 cathschlper10k secularper10k, cluster(cfips)
estimates store col5

*Ohio (using county fixed effects to show Trump vote share within school districts in the same county matter!)
oprobit openedup i.loc2 i.cfips ohiotrump distsize lninc50avgall logpps pwhite, cluster(cfips)
estimates store col6


esttab col1 col2 col3 col4 col5 col6 using table1.rtf, replace ///
	se pr2 star(† 0.1 * 0.05 ** 0.01 *** 0.001)  ///
	keep(trump distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k pac2018 cbmoe cut1 cut2) ///
	rename(ohiotrump "trump" ptrump "trump") ///
	coeflabels(trump "Partisanship (Trump vote)") ///
	order(trump meancasesper2 cathschlper10k secularper10k distsize cbmoe pac2018  cathschlper10k secularper10k lninc50avgall logpps pwhite cut1 cut2  ) ///
	indicate("State FE = *st" "Locale FE = 2.loc2" "State X Locale FE = *st#1.loc2"  "Commuter Zones FE = *czi2000" "County FE = *cfips") ///
	nonotes	///
	addnotes("† p<0.1, * p<0.05, ** p<0.01, *** p<0.001" "Clustered Standard errors in parentheses.") ///
	mgroup("National" "Washington State" "Ohio", pattern(1 0 0 0 1  1)) ///
	mlabels(none) eqlabels(none) label
	
	
	
	
************************************
************************************
************************************
************************************
************************************
*Code to reproduce Figure 2
************************************
************************************
************************************
************************************
************************************

use "Main analysis file.dta", clear

probit fullyopen i.st i.loc2 distsize logpps lninc50avgall pwhite ptrump  meancasesper2  cathschlper10k secularper10k, cluster(st) 
margins, at(ptrump=(.2(.2)1)) atmeans saving(margopen1, replace)
probit fullyremote i.st i.loc2 distsize logpps lninc50avgall pwhite ptrump  meancasesper2  cathschlper10k secularper10k, cluster(st) 
margins, at(ptrump=(.2(.2)1)) atmeans saving(margremote1, replace)
combomarginsplot margopen1 margremote1, labels("In-Person Learning" "Fully Remote/Closed") xtitle("Trump Vote Share")  ytitle("") name(m1, replace)

probit fullyopen i.st i.loc2 distsize logpps lninc50avgall pwhite ptrump  meancasesper2  cathschlper10k secularper10k, cluster(st) 
margins, at(meancasesper2=(0(10)50)) atmeans   saving(margopen2, replace)
probit fullyremote i.st i.loc2 distsize logpps lninc50avgall pwhite ptrump  meancasesper2  cathschlper10k secularper10k, cluster(st) 
margins, at(meancasesper2=(0(10)50)) atmeans   saving(margremote2, replace)
combomarginsplot margopen2 margremote2, labels("In-Person Learning" "Fully Remote/Closed") xtitle("COVID case rate") ytitle("")name(m2, replace)

probit fullyopen i.st i.loc2 distsize logpps lninc50avgall pwhite ptrump  meancasesper2  cathschlper10k secularper10k, cluster(st) 
margins, at(distsize=(0(1.5)13))  atmeans   saving(margopen3, replace)
probit fullyremote i.st i.loc2 distsize logpps lninc50avgall pwhite ptrump  meancasesper2  cathschlper10k secularper10k, cluster(st) 
margins, at(distsize=(0(1.5)13)) atmeans   saving(margremote3, replace)
combomarginsplot margopen3 margremote3, labels("In-Person Learning" "Fully Remote/Closed") xtitle("District size") ytitle("") name(m3, replace)

probit fullyopen i.st i.loc2 distsize logpps lninc50avgall pwhite ptrump  meancasesper2  cathschlper10k secularper10k, cluster(st) 
margins, at(cathschlper10k=(0(2)8)) atmeans    saving(margopen4, replace)
probit fullyremote i.st i.loc2 distsize logpps lninc50avgall pwhite ptrump  meancasesper2  cathschlper10k secularper10k, cluster(st) 
margins, at(cathschlper10k=(0(2)8)) atmeans   saving(margremote4, replace)
combomarginsplot margopen4 margremote4, labels("In-Person Learning" "Fully Remote/Closed") xtitle("Catholic schools per-capita") ytitle("") name(m4, replace)
grc1leg m1 m2 m3 m4, scale(0.6) ycommon imargin(4 4 4 4) cols(2) scheme(s1manual)
graph export "Figure 2", replace as(png)






************************************
************************************
************************************
************************************
************************************
*Code to reproduce Figure 3
************************************
************************************
************************************
************************************
************************************

use "Main analysis file.dta", clear


* Figure 3: showing marginal effects of a standard deviation increase in all public health variables

** FIGURE ***

*meancasesper1
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump  meancasesper1  cathschlper10k secularper10k, cluster(st) 
summarize meancasesper1
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(meancasesper1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store per1remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump  meancasesper1  cathschlper10k secularper10k, cluster(st) 
margins, at(meancasesper1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store per1open

*meancasesper2
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   meancasesper2 cathschlper10k secularper10k, cluster(st) 
summarize meancasesper2
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(meancasesper2 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store per2remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump  meancasesper2  cathschlper10k secularper10k, cluster(st) 
margins, at(meancasesper2 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store per2open

*meancasesper_aug_all
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   meancasesper_aug_all cathschlper10k secularper10k, cluster(st) 
summarize meancasesper_aug_all
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(meancasesper_aug_all = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store aug_allremote
 
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   meancasesper_aug_all cathschlper10k secularper10k, cluster(st) 
margins, at(meancasesper_aug_all = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store aug_allopen
 
*cumcasesper8_1
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump  cumcasesper8_1  cathschlper10k secularper10k, cluster(st) 
summarize cumcasesper8_1
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(cumcasesper8_1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store cum_8_1remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump  cumcasesper8_1  cathschlper10k secularper10k, cluster(st) 
margins, at(cumcasesper8_1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store cum_8_1open
 
*cumcasesper8_31
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   cumcasesper8_31 cathschlper10k secularper10k, cluster(st) 
summarize cumcasesper8_31
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(cumcasesper8_31 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store cum_8_31remote
 
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   cumcasesper8_31 cathschlper10k secularper10k, cluster(st) 
margins, at(cumcasesper8_31 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store cum_8_31open
 
*meandeathsper1
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump  meandeathsper1  cathschlper10k secularper10k, cluster(st) 
summarize meandeathsper1
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(meandeathsper1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store deaths1remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump  meandeathsper1  cathschlper10k secularper10k, cluster(st)  
 margins, at(meandeathsper1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store deaths1open
 
*meandeathsper2
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   meandeathsper2 cathschlper10k secularper10k, cluster(st) 
summarize meandeathsper2
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(meandeathsper2 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store deaths2remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   meandeathsper2 cathschlper10k secularper10k, cluster(st) 
margins, at(meandeathsper2 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store deaths2open

*meandeathsper_aug_all
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   meandeathsper_aug_all cathschlper10k secularper10k, cluster(st) 
summarize meandeathsper_aug_all
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(meandeathsper_aug_all = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store deaths_aug_allremote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   meandeathsper_aug_all cathschlper10k secularper10k, cluster(st) 
margins, at(meandeathsper_aug_all = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store deaths_aug_allopen

*cumdeathsper8_1
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump  cumdeathsper8_1  cathschlper10k secularper10k, cluster(st) 
summarize cumdeathsper8_1
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(cumdeathsper8_1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store cumdeaths_8_1remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump  cumdeathsper8_1  cathschlper10k secularper10k, cluster(st) 
margins, at(cumdeathsper8_1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store cumdeaths_8_1open

*cumdeathsper8_31
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   cumdeathsper8_31 cathschlper10k secularper10k, cluster(st) 
summarize cumdeathsper8_31
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(cumdeathsper8_31 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store cumdeaths_8_31remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   cumdeathsper8_31 cathschlper10k secularper10k, cluster(st) 
margins, at(cumdeathsper8_31 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store cumdeaths_8_31open

*meancombo1
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump  meancombo1  cathschlper10k secularper10k, cluster(st) 
summarize meancombo1
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(meancombo1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store combo1remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump  meancombo1  cathschlper10k secularper10k, cluster(st) 
margins, at(meancombo1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store combo1open

*meancombo2
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   meancombo2 cathschlper10k secularper10k, cluster(st) 
summarize meancombo2
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(meancombo2 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store combo2remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   meancombo2 cathschlper10k secularper10k, cluster(st) 
margins, at(meancombo2 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store combo2open

*meancombo_aug_all
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   meancombo_aug_all cathschlper10k secularper10k, cluster(st) 
summarize meancombo_aug_all
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(meancombo_aug_all = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store combo_aug_allremote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   meancombo_aug_all cathschlper10k secularper10k, cluster(st) 
margins, at(meancombo_aug_all = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store combo_aug_allopen

*hosp1
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp1 cathschlper10k secularper10k, cluster(st) 
summarize hosp1
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(hosp1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store hosp1remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp1 cathschlper10k secularper10k, cluster(st) 
margins, at(hosp1 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store hosp1open

*hosp2
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp2 cathschlper10k secularper10k, cluster(st) 
summarize hosp2
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(hosp2 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store hosp2remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp2 cathschlper10k secularper10k, cluster(st) 
margins, at(hosp2 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store hosp2open

*hosp3
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp3 cathschlper10k secularper10k, cluster(st) 
summarize hosp3
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(hosp3 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store hosp3remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp3 cathschlper10k secularper10k, cluster(st) 
margins, at(hosp3 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store hosp3open


*hosp4
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp4 cathschlper10k secularper10k, cluster(st) 
summarize hosp4
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(hosp4 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store hosp4remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp4 cathschlper10k secularper10k, cluster(st) 
margins, at(hosp4 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store hosp4open

*hosp5
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp5 cathschlper10k secularper10k, cluster(st) 
summarize hosp5
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(hosp5 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store hosp5remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp5 cathschlper10k secularper10k, cluster(st) 
margins, at(hosp5 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store hosp5open

*hosp10
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp10 cathschlper10k secularper10k, cluster(st) 
summarize hosp10
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(hosp10 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store hosp10remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp10 cathschlper10k secularper10k, cluster(st) 
margins, at(hosp10 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store hosp10open

*hosp11
quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp11 cathschlper10k secularper10k, cluster(st) 
summarize hosp11
local covid = r(mean)
local covid2 = r(mean) + r(sd)
margins, at(hosp11 = (`covid' `covid2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store hosp11remote

quietly oprobit openedup i.st i.loc2 distsize lninc50avgall logpps pwhite ptrump   hosp11 cathschlper10k secularper10k, cluster(st) 
margins, at(hosp11 = (`covid' `covid2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store hosp11open

*Baseline Nationwide estimation (using state and locale FEs)
quietly oprobit openedup i.st i.loc2 ptrump distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k, cluster(st) 
summarize ptrump
local trump = r(mean)
local trump2 = r(mean) + r(sd)
margins, at(ptrump = (`trump' `trump2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store ptrumpremote

quietly oprobit openedup i.st i.loc2 ptrump distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k, cluster(st) 
margins, at(ptrump = (`trump' `trump2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store ptrumpopen

*Baseline Nationwide estimation (using state and locale FEs)
quietly oprobit openedup i.st i.loc2 ptrump distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k, cluster(st) 
summarize cathschlper10k
local cath = r(mean)
local cath2 = r(mean) + r(sd)
margins, at(cathschlper10k = (`cath' `cath2') ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store cathremote

quietly oprobit openedup i.st i.loc2 ptrump distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k, cluster(st) 
margins, at(cathschlper10k = (`cath' `cath2') ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store cathopen

*Collective Bargaining
oprobit openedup i.st i.loc2 ptrump cbmoe distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k, cluster(st)
margins, at(cbmoe = (0 1) ) post coeflegend
lincomest  _b[1bn._predict#2._at]-_b[1bn._predict#1bn._at]
est store cbremote

oprobit openedup i.st i.loc2 ptrump cbmoe distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k, cluster(st)
margins, at(cbmoe = (0 1) ) post coeflegend
lincomest  _b[3._predict#2._at]-_b[3._predict#1bn._at]
est store cbopen

coefplot (per1open,  mcol(gray) rename((1)="Case Rates Early Aug.") label("Fully Open") msymbol(smcircle) )(per1remote, mcol(gray) msymbol(X)  rename((1)="Case Rates Early Aug.") label("Fully Remote")) ///
(per2remote, mcol(gray) msymbol(X)  rename((1)="Case Rates Late Aug.") nokey) (per2open,  msymbol(smcircle) rename((1)="Case Rates Late Aug.") mcol(gray) nokey) ///
(aug_allremote,  mcol(gray) msymbol(X)  rename((1)="Case Rates Aug. All") nokey) (aug_allopen,  msymbol(smcircle) rename((1)="Case Rates Aug. All")  mcol(gray) nokey) ///
(deaths1remote,  mcol(gray) msymbol(X)  rename((1)="Death Rates Early Aug.") nokey) (deaths1open,  msymbol(smcircle) rename((1)="Death Rates Early Aug.") mcol(gray) nokey) ///
(deaths2remote,  mcol(gray) msymbol(X)  rename((1)="Death Rates Late Aug.") nokey) (deaths2open,  msymbol(smcircle) rename((1)="Death Rates Late Aug.") mcol(gray) nokey) ///
(deaths_aug_allremote,  mcol(gray) msymbol(X)  rename((1)="Death Rates Aug. All") nokey) (deaths_aug_allopen,  msymbol(smcircle)  rename((1)="Death Rates Aug. All") mcol(gray) nokey) ///
(cum_8_1remote,  mcol(gray) msymbol(X)  rename((1)="Total Cases on Aug. 1") nokey) (cum_8_1open,  msymbol(smcircle) rename((1)="Total Cases on Aug. 1") mcol(gray) nokey) ///
(cum_8_31remote,  mcol(gray) msymbol(X)  rename((1)="Total Cases on Aug. 31") nokey) (cum_8_31open,  msymbol(smcircle) rename((1)="Total Cases on Aug. 31") mcol(gray) nokey) ///
(cumdeaths_8_1remote,  mcol(gray) msymbol(X)  rename((1)="Total Deaths on Aug. 1") nokey) (cumdeaths_8_1open,  msymbol(smcircle) rename((1)="Total Deaths on Aug. 1") mcol(gray) nokey) ///
(cumdeaths_8_31remote,  mcol(gray) msymbol(X)  rename((1)="Total Deaths on Aug. 31") nokey) (cumdeaths_8_31open,  msymbol(smcircle) rename((1)="Total Deaths on Aug. 31") mcol(gray) nokey) ///
(combo1remote,  mcol(gray) msymbol(X)  rename((1)="Combo Indicators Early Aug.") nokey) (combo1open,  msymbol(smcircle) rename((1)="Combo Indicators Early Aug.") mcol(gray) nokey) ///
(combo2remote,  mcol(gray) msymbol(X)  rename((1)="Combo Indicators Late Aug.") nokey) (combo2open,  msymbol(smcircle) rename((1)="Combo Indicators Late Aug.") mcol(gray) nokey) ///
(combo_aug_allremote,  mcol(gray) msymbol(X)  rename((1)="Combo Indicators Aug. All") nokey) (combo_aug_allopen,  msymbol(smcircle) rename((1)="Combo Indicators Aug. All") mcol(gray) nokey) ///
(hosp1remote,  mcol(gray) msymbol(X)  rename((1)="Hospitalized 7-31-20") nokey) (hosp1open,  msymbol(smcircle) rename((1)="Hospitalized 7-31-20") mcol(gray) nokey) ///
(hosp2remote,  mcol(gray) msymbol(X)  rename((1)="Hospitalized 8-07-20") nokey) (hosp2open,  msymbol(smcircle) rename((1)="Hospitalized 8-07-20") mcol(gray) nokey) ///
(hosp3remote,  mcol(gray) msymbol(X)  rename((1)="Hospitalized 8-14-20") nokey) (hosp3open,  msymbol(smcircle) rename((1)="Hospitalized 8-14-20") mcol(gray) nokey) ///
(hosp4remote,  mcol(gray) msymbol(X)  rename((1)="Hospitalized 8-21-20") nokey) (hosp4open,  msymbol(smcircle) rename((1)="Hospitalized 8-21-20") mcol(gray) nokey) ///
(hosp5remote,  mcol(gray) msymbol(X)  rename((1)="Hospitalized 8-28-20") nokey) (hosp5open,  msymbol(smcircle) rename((1)="Hospitalized 8-28-20") mcol(gray) nokey) ///
(hosp10remote,  mcol(gray) msymbol(X)  rename((1)="Hospitalized 10-2-20") nokey) (hosp10open,  msymbol(smcircle) rename((1)="Hospitalized 10-2-20") mcol(gray) nokey) ///
(hosp11remote,  mcol(gray) msymbol(X)  rename((1)="Hospitalized 10-9-20") nokey) (hosp11open,  msymbol(smcircle) rename((1)="Hospitalized 10-9-20") mcol(gray) nokey) ///
(ptrumpremote,  mcol(gray) msymbol(X)  rename((1)="Partisanship (Trump vote)") nokey) (ptrumpopen,  msymbol(smcircle) rename((1)="Partisanship (Trump vote)") mcol(gray) nokey) ///
(cbremote,  mcol(gray) msymbol(X)  rename((1)="Collective bargaining") nokey) (cbopen,  msymbol(smcircle) rename((1)="Collective bargaining") mcol(gray) nokey) ///
(cathremote,  mcol(gray) msymbol(X)  rename((1)="Catholic Schools") nokey) (cathopen,  msymbol(smcircle) rename((1)="Catholic Schools") mcol(gray) nokey) ///
, ciopts(lcolor(gray)) xline(0) ysize(7) xsize(10)
graph export "Figure 3", replace as(png)





************************************
************************************
************************************
************************************
************************************
*Code to replicate Figure 4
************************************
************************************
************************************
************************************
************************************

use "Main analysis file.dta", clear

probit fullyopen i.st i.loc2 distsize  lninc50avgall pwhite c.ptrump##c.meancasesper2 cathschlper10k secularper10k, cluster(st)
margins, at(meancasesper2=(0 10 20 30 40 50) ptrump=(.25 .5 .75)) vsquish
marginsplot, name(open, replace) 

probit fullyremote i.st i.loc2 distsize  lninc50avgall pwhite c.ptrump##c.meancasesper2  cathschlper10k secularper10k, cluster(st)
margins, at(meancasesper2=(0 10 20 30 40 50) ptrump=(.25 .5 .75)) vsquish
marginsplot, name(closed, replace) 
grc1leg open closed, ycommon scale(0.8) imargin(4 4 4 4) cols(2) scheme(s1manual) 
graph export "Figure 4", replace as(png)







************************************
************************************
************************************
************************************
************************************
*Code to replicate Table A1
************************************
************************************
************************************
************************************
************************************

use "Main analysis file.dta", clear

labsumm ptrump ohiotrump distsize lninc50avgall logpps ///
pwhite cathschlper10k secularper10k pac2018 cbmoe meancasesper1 ///
meancasesper2 meancasesper_aug_all meandeathsper1 meandeathsper2 ///
meandeathsper_aug_all cumcasesper8_1 cumcasesper8_31 cumdeathsper8_1 ///
cumdeathsper8_31 meancombo1 meancombo2 meancombo_aug_all hosp1-hosp5 hosp10 hosp11


************************************
************************************
************************************
************************************
************************************
*Code to replicate Table A2
************************************
************************************
************************************
************************************
************************************

use "Main analysis file.dta", clear

oprobit openedup i.st i.loc2 c.ptrump##c.meancasesper2 distsize lninc50avgall logpps pwhite cathschlper10k secularper10k, cluster(st)
estimates store a3

esttab a3 using tablea2.rtf, ///
	se pr2 star(† 0.1 * 0.05 ** 0.01 *** 0.001)  ///
	keep(ptrump c.ptrump#c.meancasesper2 distsize lninc50avgall logpps pwhite meancasesper2 cathschlper10k secularper10k cut1 cut2) ///
	order(ptrump  meancasesper2 c.ptrump#c.meancasesper2 cathschlper10k secularper10k distsize cathschlper10k secularper10k lninc50avgall logpps pwhite cut1 cut2  ) ///
	coeflabels(c.ptrump#c.meancasesper2 "Partisanship (Trump vote) X COVID case rate") ///
	indicate("State FE = *st" "Locale FE = *loc2") ///
	nonotes	///
	addnotes("† p<0.1, * p<0.05, ** p<0.01, *** p<0.001" "Clustered Standard errors in parentheses.") ///
	mlabels(none) eqlabels(none)  label


erase margopen1.dta 
erase margopen2.dta
erase margopen3.dta 
erase margopen4.dta 
erase margremote1.dta 
erase margremote2.dta 
erase margremote3.dta 
erase margremote4.dta


